#include <feel/options.hpp>
#include <feel/feelalg/backend.hpp>
#include <feel/feeldiscr/functionspace.hpp>
#include <feel/feeldiscr/region.hpp>
#include <feel/feelpoly/im.hpp>
#include <feel/feelfilters/gmsh.hpp>
#include <feel/feelfilters/exporter.hpp>
#include <feel/feelpoly/polynomialset.hpp>
#include <feel/feelvf/vf.hpp>

#include <feel/feelfilters/geotool.hpp>

//#include "../Schwarz2/loadmesh.hpp"

/** use Feel namespace */
using namespace Feel;
using namespace Feel::vf;

//---------------------------------------------------------------------//
 

# define INTERP_OPERATOR						\
  BOOST_PP_TUPLE_TO_LIST( 2, ( ( idv   , "idv"  , *,"*", idv, "idv",1 ), \
			       ( gradv , "gradv", *,"*", gradv , "gradv"  ,0 ) ) \
			  )						\
  /**/

// Accessors for the operator datatype.
# define INTERP_OPERATOR_SYMBOL(O) BOOST_PP_TUPLE_ELEM(7, 0, O)
# define INTERP_OPERATOR_NAME(O) BOOST_PP_TUPLE_ELEM(7, 1, O)
# define INTERP_OPERATOR_OPERATION(O) BOOST_PP_TUPLE_ELEM(7, 2, O)
# define INTERP_OPERATOR_OP_STR(O) BOOST_PP_TUPLE_ELEM(7, 3, O)
# define INTERP_OPERATOR_SYMBOL2(O) BOOST_PP_TUPLE_ELEM(7, 4, O)
# define INTERP_OPERATOR_STR2(O) BOOST_PP_TUPLE_ELEM(7, 5, O)
# define INTERP_OPERATOR_ID(O) BOOST_PP_TUPLE_ELEM(7, 6, O)


//---------------------------------------------------------------------//

#define INTERP_CHAMP					\
  BOOST_PP_TUPLE_TO_LIST( 3, ( ( u1 , "u1",mesh1 ),		\
			       ( u2 , "u2",mesh2 ),			\
			       ( u_2p , "u_2p",meshBase_2p ) )			\
			  )				\
  /**/

# define INTERP_CHAMP_SYMBOL(O) BOOST_PP_TUPLE_ELEM(3, 0, O)
# define INTERP_CHAMP_NAME(O) BOOST_PP_TUPLE_ELEM(3, 1, O)
# define INTERP_CHAMP_MESH(O) BOOST_PP_TUPLE_ELEM(3, 2, O)


//---------------------------------------------------------------------//
#define INTERP_MESH_PARTION ( meshBase_2p, mesh1, mesh2 )

#define INTERP_MESH							\
  BOOST_PP_TUPLE_TO_LIST( 3, ( ( mesh1       , "mesh1"       , u1  , "u1"  , 1 ), \
			       ( mesh2       , "mesh2"       , u2  , "u2"  , 2 ), \
			       ( meshBase_2p , "meshBase_2p" , u_2p, "u_2p", 0 ) ) \
			  )						\
  /**/

# define INTERP_MESH_SYMBOL(O) BOOST_PP_TUPLE_ELEM(5, 0, O)
# define INTERP_MESH_NAME(O) BOOST_PP_TUPLE_ELEM(5, 1, O)
# define INTERP_MESH_CHAMPINMESH2(O) BOOST_PP_TUPLE_ELEM(5, 2, O)
# define INTERP_MESH_CHAMPINMESH(O) BOOST_PP_TUPLE_ELEM(5, 3, O)
# define INTERP_MESH_TYPEOFPARTITION(O) BOOST_PP_TUPLE_ELEM(5, 4, O)

//---------------------------------------------------------------------//

# define INTERP_TEST_EXECUTION_ELEMENTSMESH(_, MOT)	\
  INTERP_TEST_EXECUTION_ELEMENTSMESH_CODE MOT		\
  /**/



#define INTERP_TEST_EXPRESSION(M,O,T)					\
  INTERP_OPERATOR_SYMBOL(O)(INTERP_CHAMP_SYMBOL(T))			\
  *									\
  BOOST_PP_IF(INTERP_OPERATOR_ID(O),					\
	      INTERP_OPERATOR_SYMBOL2(O)(INTERP_CHAMP_SYMBOL(T)) ,	\
	      trans( INTERP_OPERATOR_SYMBOL2(O)(INTERP_CHAMP_SYMBOL(T)) ) ) \
				

#define INTERP_TEST_EXECUTION_ELEMENTSMESH_CODE(M,O,T)			\
  std::cout								\
  << "\nintegrate(elements(" << INTERP_MESH_NAME(M) << ") , "		\
  << INTERP_OPERATOR_NAME(O)<<"(" << INTERP_CHAMP_NAME(T) << ")"	\
  << INTERP_OPERATOR_OP_STR(O)						\
  << BOOST_PP_IF(INTERP_OPERATOR_ID(O),					\
		 INTERP_OPERATOR_STR2(O) <<"(" << INTERP_CHAMP_NAME(T) << ")" , \
		 "trans(" << INTERP_OPERATOR_STR2(O) <<"("<<INTERP_CHAMP_NAME(T)<<")" <<")" \
		 )							\
  << ")";								\
  __time.restart();							\
  __value = std::sqrt( integrate( elements( INTERP_MESH_SYMBOL(M) ),	\
		       INTERP_TEST_EXPRESSION(M,O,T)			\
				  ).evaluate()(0,0) );			\
  __t = __time.elapsed();						\
  __MapMemoryElements[INTERP_MESH_NAME(M)][INTERP_CHAMP_NAME(T)][INTERP_OPERATOR_NAME(O)] = __value; \
  std::cout << " : " << __value << " in " << __t << "sec";		\
  /**/				


# define INTERP_TEST_EXECUTION_BOUNDARYINTERFACE(_, MOT)	\
  INTERP_TEST_EXECUTION_BOUNDARYINTERFACE_CODE MOT		\
  /**/

#define INTERP_TEST_EXECUTION_BOUNDARYINTERFACE_CODE(M,O,T)		\
  std::cout << "\nintegrate( markedfaces(" << INTERP_MESH_NAME(M) << "->markerName(Interface))," \
  << INTERP_OPERATOR_NAME(O)<<"(" << INTERP_CHAMP_NAME(T) << ")"	\
  << "*" << BOOST_PP_IF(INTERP_OPERATOR_ID(O),				\
			INTERP_OPERATOR_STR2(O) << "(" << INTERP_CHAMP_NAME(T) << ")" , \
			"trans (" << INTERP_OPERATOR_STR2(O) << "(" INTERP_CHAMP_NAME(T) << ")" ); \
  __time.restart();							\
  __value = std ::sqrt( integrate( markedfaces( INTERP_MESH_SYMBOL(M), INTERP_MESH_SYMBOL(M)->markerName("Interface") ), \
				   INTERP_TEST_EXPRESSION(M,O,T)	\
				   ).evaluate()(0,0) );			\
  __t = __time.elapsed();						\
  __MapMemoryInterface[INTERP_MESH_NAME(M)][INTERP_CHAMP_NAME(T)][INTERP_OPERATOR_NAME(O)] = __value; \
  std::cout << " : " <<  __value << " in " << __t << "sec";		\
  /**/				


//---------------------------------------------------------------------//

# define INTERP_TEST_VALIDATION_ELEMENTSMESH(_, MO)	\
  INTERP_TEST_VALIDATION_ELEMENTSMESH_CODE MO		\
   /**/


# define INTERP_TEST_VALIDATION_ELEMENTSMESH_CODE( M,O)			\
  std::cout<<"\n"<< INTERP_MESH_NAME(M) << " " << INTERP_OPERATOR_NAME(O); \
  BOOST_PP_IF( BOOST_PP_EQUAL( INTERP_MESH_TYPEOFPARTITION(M),1 ),	\
	       if ( ( std::abs( __MapMemoryElements[INTERP_MESH_NAME(M)]["u_2p"][INTERP_OPERATOR_NAME(O)] \
				- __MapMemoryElements[INTERP_MESH_NAME(M)]["u1"][INTERP_OPERATOR_NAME(O)])<1e-6) \
		    && ( __MapMemoryElements[INTERP_MESH_NAME(M)]["u2"][INTERP_OPERATOR_NAME(O)] == 0 ) ) \
		 std::cout << " TRUE"; else std::cout << " FALSE";\
	       , )						   \
  BOOST_PP_IF( BOOST_PP_EQUAL( INTERP_MESH_TYPEOFPARTITION(M),2 ),	\
	       if ( ( std::abs( __MapMemoryElements[INTERP_MESH_NAME(M)]["u_2p"][INTERP_OPERATOR_NAME(O)] \
				- __MapMemoryElements[INTERP_MESH_NAME(M)]["u2"][INTERP_OPERATOR_NAME(O)])<1e-6) \
		    && ( __MapMemoryElements[INTERP_MESH_NAME(M)]["u1"][INTERP_OPERATOR_NAME(O)] == 0 ) ) \
		 std::cout << " TRUE"; else std::cout << " FALSE";\
	       , )						   \
  BOOST_PP_IF( BOOST_PP_EQUAL( INTERP_MESH_TYPEOFPARTITION(M),0 ), \
	       if ( ( std::abs( __MapMemoryElements[INTERP_MESH_NAME(M)]["u1"][INTERP_OPERATOR_NAME(O)] \
				- __MapMemoryElements["mesh1"]["u1"][INTERP_OPERATOR_NAME(O)])<1e-6) \
		    && ( std::abs( __MapMemoryElements[INTERP_MESH_NAME(M)]["u2"][INTERP_OPERATOR_NAME(O)] \
				- __MapMemoryElements["mesh2"]["u2"][INTERP_OPERATOR_NAME(O)])<1e-6 ) ) \
		 std::cout << " TRUE"; else std::cout << " FALSE"; 	\
	       ,)							\
  /*  std::cout << " " <<  std::abs(__MapMemoryElements[INTERP_MESH_NAME(M)]["u_2p"][INTERP_OPERATOR_NAME(O)] \
  - __MapMemoryElements[INTERP_MESH_NAME(M)]["u1"][INTERP_OPERATOR_NAME(O)] \
  - __MapMemoryElements[INTERP_MESH_NAME(M)]["u2"][INTERP_OPERATOR_NAME(O)]\
  + __MapMemoryInterface[INTERP_MESH_NAME(M)][INTERP_MESH_CHAMPINMESH(M)][INTERP_OPERATOR_NAME(O)])*/ \
 ; \
  /**/

# define INTERP_TEST_VALIDATION_BOUNDARYINTERFACE(_, MOT)	\
  INTERP_TEST_VALIDATION_BOUNDARYINTERFACE_CODE MOT		\
   /**/


# define INTERP_TEST_VALIDATION_BOUNDARYINTERFACE_CODE( M,O,T)				\
  std::cout<<"\n"<< INTERP_MESH_NAME(M) << " " << INTERP_OPERATOR_NAME(O) << " " << INTERP_CHAMP_NAME(T); \
  if ( std::abs( __MapMemoryInterface[INTERP_MESH_NAME(M)][INTERP_CHAMP_NAME(T)][INTERP_OPERATOR_NAME(O)] \
		- __MapMemoryInterface[INTERP_MESH_NAME(M)][INTERP_MESH_CHAMPINMESH(M)][INTERP_OPERATOR_NAME(O)]) < 1e-6 ) \
    std::cout << " TRUE";						\
  else std::cout << " FALSE";						\
  std::cout<< "  " << ( std::abs(__MapMemoryInterface[INTERP_MESH_NAME(M)][INTERP_CHAMP_NAME(T)][INTERP_OPERATOR_NAME(O)] \
				 - __MapMemoryInterface[INTERP_MESH_NAME(M)][INTERP_MESH_CHAMPINMESH(M)][INTERP_OPERATOR_NAME(O)]) ); \
  /**/


	

inline
po::options_description
makeOptions()
{
  po::options_description code1doptions("Code1D options");
  code1doptions.add_options()
    ("hsize1", po::value<double>()->default_value( 0.1 ), "mesh size1")
    ("hsize2", po::value<double>()->default_value( 0.1 ), "mesh size2")
    ("gmsh", po::value<float>()->default_value( 2.1 ), " version of gmsh( 2.0 or 2.1)")
    ;
  return code1doptions.add( Feel::feel_options() );
}

/**
 * This routine defines some information about the application like
 * authors, version, or name of the application. The data returned is
 * typically used as an argument of a Feel::Application subclass.
 *
 * \return some data about the application.
 */
inline
AboutData
makeAbout()
{
  AboutData about( "Test_Interpolate" ,
		   "Test_Interpolate" ,
		   "0.2",
		   "Code 1D for -Delta u = f",
		   Feel::AboutData::License_GPL,
		   "Copyright (c) 2010 Universite Joseph Fourier");

  about.addAuthor("Vincent Chabannes", "developer", "vincent.chabannes@imag.fr", "");
  return about;

}


template<uint Dim,uint N,uint OrderGeo>
class Localisation
  :
  public Application
{
  typedef Application super;
public:

  //! Polynomial order \f$P_2\f$
  //    static const uint16_type Order = 5;
  ///static const uint Order = N;

  typedef double value_type;

  typedef Backend<value_type> backend_type;
  typedef boost::shared_ptr<backend_type> backend_ptrtype;

  typedef backend_type::sparse_matrix_type sparse_matrix_type;
  typedef backend_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
  typedef backend_type::vector_type vector_type;
  typedef backend_type::vector_ptrtype vector_ptrtype;

  typedef Simplex<Dim,OrderGeo,Dim> convex_type;
  typedef Mesh<convex_type> mesh_type;
  typedef boost::shared_ptr<mesh_type> mesh_ptrtype;

  typedef bases<Lagrange<N,Scalar,PointSetFekete> > basis_type;
  //  typedef bases<Lagrange<N,Scalar/*,PointSetFekete*/> > basis_type;
  typedef FunctionSpace<mesh_type, basis_type> space_type;
  typedef boost::shared_ptr<space_type> space_ptrtype;
  typedef typename space_type::element_type element_type;

  typedef Exporter<mesh_type> export_type;
  typedef boost::shared_ptr<export_type> export_ptrtype;

  /**
   * Constructor
   */
  Localisation( int argc, char** argv, AboutData const& ad,
		po::options_description const& od )
    :
    super( argc, argv, ad, od ),
    M_backend( backend_type::build( this->vm() ) ),
    //M_backend( backend_type::build( BACKEND_TRILINOS ) ),
    meshSize1( this->vm()["hsize1"].template as<double>() ),
    meshSize2( this->vm()["hsize2"].template as<double>() ),
    exporter( Exporter<mesh_type>::New( this->vm(), this->about().appName() ) )
  {
  }

  mesh_ptrtype createMesh( std::string,std::string );

  mesh_ptrtype createMesh( double meshSize );

  mesh_ptrtype createMesh2( double meshSize );
    

  void run();


private:

  void testEvaluate(boost::tuple<mesh_ptrtype,mesh_ptrtype> __mesh);
  void testAssemblage(boost::tuple<mesh_ptrtype,mesh_ptrtype> __mesh);

  void exportResults( element_type& u,std::string __str = "u" );

private:

  // linear algebra backend
  backend_ptrtype M_backend;

  // mesh characteristic size
  double meshSize1;
  double meshSize2;

  // exporter factory
  export_ptrtype exporter;
}; // Interpolate



template<uint Dim,uint N,uint OrderGeo>
typename Localisation<Dim,N,OrderGeo>::mesh_ptrtype
Localisation<Dim,N,OrderGeo>::createMesh( std::string __name,std::string __str) 
{
  mesh_ptrtype mesh( new mesh_type );

  Gmsh gmsh;
  gmsh.setOrder(OrderGeo);

  std::string fname = gmsh.generate( __name, __str,false,false,false  );

  ImporterGmsh<mesh_type> import( fname );
    
  float gmsh_version = this->vm()["gmsh"].template as<float>();
  if (gmsh_version==2.0) import.setVersion( "2.0" );
  else if (gmsh_version==2.1) import.setVersion( "2.1" );

  mesh->accept( import );

  return mesh;

}


template<uint Dim,uint N,uint OrderGeo>
typename Localisation<Dim,N,OrderGeo>::mesh_ptrtype
Localisation<Dim,N,OrderGeo>::createMesh( double meshSize ) {
  mesh_ptrtype mesh( new mesh_type );

  std::ostringstream ostr;
    
  GeoTool::Node x1(-1,-1);
  GeoTool::Node x2(0.2,1);
  GeoTool::Rectangle R1(0.1,"LEFT",x1,x2);
  R1.setMarker(_type="line",_name="Interface",_marker2=true);
  R1.setMarker(_type="line",_name="Bord",_marker1=true,_marker3=true,_marker4=true);
  R1.setMarker(_type="surface",_name="Omega1",_marker1=true);
   
  Gmsh gmsh;

  gmsh.setOrder(OrderGeo);

  std::string fname = gmsh.generate( "domain", R1.geoStr(),false,false,false  );

  ImporterGmsh<mesh_type> import( fname );
    
  float gmsh_version = this->vm()["gmsh"].template as<float>();

  if (gmsh_version==2.0) import.setVersion( "2.0" );
  else if (gmsh_version==2.1) import.setVersion( "2.1" );

  mesh->accept( import );

  return mesh;

}



template<uint Dim,uint N,uint OrderGeo>
void
Localisation<Dim,N,OrderGeo>::exportResults( element_type& U, std::string __str )
{

  export_ptrtype exporterLoc( Exporter<mesh_type>::New( this->vm(), this->about().appName()+__str ) );

  Log() << "exportResults starts\n";
  exporterLoc->step(0)->setMesh( U.functionSpace()->mesh() );
  exporterLoc->step(0)->add( __str, U );
  exporterLoc->save();
  Log() << "exportResults finish\n";
} // Interpolate::export


template<uint Dim,uint N,uint OrderGeo>
void
Localisation<Dim,N,OrderGeo>::testEvaluate(boost::tuple<mesh_ptrtype,mesh_ptrtype> __mesh )
{

  mesh_ptrtype mesh1 = boost::get<0>(__mesh);
  mesh_ptrtype mesh2 = boost::get<1>(__mesh);

  space_ptrtype Xh1 = space_type::New( mesh1 );
  space_ptrtype Xh2 = space_type::New( mesh2 );

  element_type u1( Xh1, "u1" );
  element_type u2( Xh2, "u2" );


  AUTO (e ,exp(Px()*Py())*sin(2*M_PI*Px()));
  AUTO (f , sin(2*M_PI*Px())*sin(2*M_PI*Py()));
  //AUTO (g , sin(0.5*M_PI*Px())*(cos(0.5*M_PI*Py())));
  //AUTO (g , Px()*(Px()-1)*Py()*(Py()-1) );

  //AUTO (h , Px()*(Py()+1) );
  AUTO (h , (Px()+1)*(Px()-1)*(Py()+1)*(Py()-1) );
  //AUTO (h , cst(1.)*Px());

  AUTO (g , 2.5*Px()+0.5);

  u1 = project( Xh1, elements(mesh1), h );
  u2 = project( Xh2, elements(mesh2), h );

  //-----------------------------------------------------------------------------------//
  /*
  double  __eval_m1_id11 = integrate( markedfaces(mesh1, mesh1->markerName("Interface")), idv(u1) ).evaluate()(0,0);

  double  __eval_m2_id11 = integrate( markedfaces(mesh2, mesh2->markerName("Interface")), idv(u1) ).evaluate()(0,0);
  //double __errId = (__eval_m1_id11 - __eval_m2_id11);///__eval_m1_id11;
double __errId = std::sqrt( std::pow(__eval_m1_id11 - __eval_m2_id11,2));///__eval_m1_id11;


 double  __eval_m1_grad11 = integrate(markedfaces(mesh1, mesh1->markerName("Interface")),
						     gradv(u1) ).evaluate()(0,0);
 double  __eval_m2_grad11 = integrate(markedfaces(mesh2, mesh2->markerName("Interface")), 
						     gradv(u1) ).evaluate()(0,0);
 double __errGrad = std::sqrt( std::pow(__eval_m1_grad11 - __eval_m2_grad11,2));//__eval_m1_grad11;
 
 double  __eval_m1_id2 = integrate( markedfaces(mesh1, mesh1->markerName("Interface")), idv(u2) ).evaluate()(0,0);
 double  __eval_m2_id2 = integrate( markedfaces(mesh2, mesh2->markerName("Interface")), idv(u2) ).evaluate()(0,0);
 double __errId2 = std::sqrt(std::pow(__eval_m1_id2 - __eval_m2_id2,2));//__eval_m1_id2;

   double  __eval_m1_grad2 = std::sqrt(integrate( markedfaces(mesh1, mesh1->markerName("Interface")),
						 gradv(u2)*gradv(u2) ).evaluate()(0,0));
  double  __eval_m2_grad2 = std::sqrt(integrate( markedfaces(mesh2, mesh2->markerName("Interface")),
						 gradv(u2)*gradv(u2) ).evaluate()(0,0));
  double __errGrad2 = (__eval_m1_grad2 - __eval_m2_grad2)/__eval_m1_grad2;

  */


  double  __errId = std::sqrt(integrate( markedfaces(mesh1, mesh1->markerName("Interface")),
						 (idv(u1)-idv(u2))*(idv(u1)-idv(u2)) ).evaluate()(0,0));

  double  __errGrad = std::sqrt(integrate( markedfaces(mesh1, mesh1->markerName("Interface")),
						 (gradv(u1)-gradv(u2))*trans(gradv(u1)-gradv(u2)) ).evaluate()(0,0));

  double  __errId2 = std::sqrt(integrate( markedfaces(mesh2, mesh2->markerName("Interface")),
						 (idv(u1)-idv(u2))*(idv(u1)-idv(u2)) ).evaluate()(0,0));

  double  __errGrad2 = std::sqrt(integrate( markedfaces(mesh2, mesh2->markerName("Interface")),
						 (gradv(u1)-gradv(u2))*trans(gradv(u1)-gradv(u2)) ).evaluate()(0,0));


  std::cout << "\n[testBoundary] idv(u1) error : " << __errId;
  std::cout << "\n[testBoundary] gradv(u1) error : " << __errGrad;
  std::cout << "\n[testBoundary] idv(u2) error : " << __errId2;
  std::cout << "\n[testBoundary] gradv(u2) error : " << __errGrad2;

}

template<uint Dim,uint OrderChamp,uint OrderGeo>
void
Localisation<Dim,OrderChamp,OrderGeo>::run()
{

  if ( this->vm().count( "help" ) )
    {
      std::cout << this->optionsDescription() << "\n";
      return;
    }

  using namespace Feel::vf;  
  //-----------------------------------------------------------------------------------//

  int MyOrder=OrderChamp;

  this->changeRepository( boost::format( "%1%/%2%/P%3%/h_%4%/" )
			  % this->about().appName()
			  % convex_type::name()
			  % MyOrder
			  % this->vm()["hsize1"].template as<double>()
			  );

  //-----------------------------------------------------------------------------------//

  //Geometry
  GeoTool::Node x1(-1,-1);
  GeoTool::Node x2(0,1);
  GeoTool::Rectangle R1( meshSize1,"LEFT",x1,x2);
  R1.setMarker(_type="line",_name="Interface",_marker2=true);
  R1.setMarker(_type="line",_name="Bord",_marker1=true,_marker3=true,_marker4=true);
  R1.setMarker(_type="surface",_name="Omega1",_marker1=true);
 
  GeoTool::Node x3(0,-1);
  GeoTool::Node x4(1,1);
  GeoTool::Rectangle R2(meshSize2,"RIGHT",x3,x4);
  R2.setMarker(_type="line",_name="Interface",_marker4=true);
  R2.setMarker(_type="line",_name="Bord",_marker1=true,_marker2=true,_marker3=true);
  R2.setMarker(_type="surface",_name="Omega2",_marker1=true);

  //-----------------------------------------------------------------------------------//

  GeoTool::Node x21(0,0);
  GeoTool::Node x22(1,1);
  GeoTool::Rectangle R(meshSize1,"RIGHT",x21,x22);
  R.setMarker(_type="line",_name="Bord",_marker1=true,_marker2=true,_marker3=true,_marker4=true);
  R.setMarker(_type="surface",_name="Omega3",_marker1=true);
  
  GeoTool::Node x_centre(0.5,0.5);
  GeoTool::Node x_bord(0.7,0.5);
  GeoTool::Circle C(meshSize2,"uncercle",x_centre,x_bord);
  C.setMarker(_type="line",_name="Interface",_marker1=true);
  //C.setMarker(_type="surface",_name="Omega2",_marker1=true);

  GeoTool::Node x2_centre(0.5,0.5);
  GeoTool::Node x2_bord(0.7,0.5);
  GeoTool::Circle C2(meshSize2,"uncercle",x2_centre,x2_bord);
  C2.setMarker(_type="line",_name="Interface",_marker1=true);
  C2.setMarker(_type="surface",_name="Omega4",_marker1=true);

  //std::cout<<(R-C).geoStr();

  //-----------------------------------------------------------------------------------//

  //Mesh
  mesh_ptrtype mesh1 = createMesh( "mesh1",R1.geoStr() );
  mesh_ptrtype mesh2 = createMesh( "mesh2",R2.geoStr() );

  mesh_ptrtype mesh3 = createMesh( "mesh3",(R-C).geoStr() );
  mesh_ptrtype mesh4 = createMesh( "mesh4",C2.geoStr() );

  std::list<boost::tuple<mesh_ptrtype,mesh_ptrtype> > __listMesh;
  __listMesh.push_back(boost::make_tuple(mesh1,mesh2));
  __listMesh.push_back(boost::make_tuple(mesh3,mesh4));

  typename std::list<boost::tuple<mesh_ptrtype,mesh_ptrtype> >::iterator itMesh = __listMesh.begin();
  typename std::list<boost::tuple<mesh_ptrtype,mesh_ptrtype> >::iterator itMesh_end = __listMesh.end();

  int i=1;
  for ( ; itMesh!=itMesh_end ; ++itMesh, ++i)
    {
      std::cout << "\n[run] start test interpolation mesh " << i;
      testEvaluate( *itMesh);
      std::cout << "\n[run] finish test interpolation mesh " << i;
    }


  /*
  mesh_ptrtype mesh1bis = createMesh( "mesh1bis",R1.geoStr() );
  space_ptrtype Xh1 = space_type::New( mesh1 );
  space_ptrtype Xh1bis = space_type::New( mesh1bis );

  element_type u1( Xh1, "u1" );
  element_type u1bis( Xh1bis, "u1bis" );

  AUTO (h , (Px()+1)*(Px()-1)*(Py()+1)*(Py()-1) );
  //AUTO (h ,exp(Px()*Py())*sin(2*M_PI*Px()));
  //AUTO (h , cst(1.) );


  u1 = project( Xh1, elements(mesh1), h );
  u1bis = project( Xh1bis, elements(mesh1bis), h );

  double __errId = 
    std::sqrt(integrate( elements(mesh1), idv(u1)*idv(u1) ).evaluate()(0,0))
    -std::sqrt(integrate( elements(mesh1bis), idv(u1)*idv(u1) ).evaluate()(0,0));

  double __errGrad = 
    std::sqrt(integrate( elements(mesh1), gradv(u1)*trans(gradv(u1)) ).evaluate()(0,0))
    -std::sqrt(integrate( elements(mesh1bis), gradv(u1)*trans(gradv(u1)) ).evaluate()(0,0));

  std::cout << "\n -----u1 id boundary Interface :" 
	    << __errId << "\n";

  std::cout << "\n -----u1 grad boundary Interface :" 
	    << __errGrad << "\n";

  double __t;
  boost::timer __time;
  double __value;
  */
  /* 
  exportResults(u_1p,"u_1p");
  exportResults(u_2p,"u_2p");
  exportResults(u1,"u1");
  exportResults(u2,"u2");
  */


  //  exportResults(u_2p);


}



/**
 * main function: entry point of the program
 */
int
main( int argc, char** argv )
{
  //dim, N, OrderGeo
  //  Localisation<2,6,1> moncode1( argc, argv, makeAbout(), makeOptions() );

  Localisation<2,6,2> moncode2( argc, argv, makeAbout(), makeOptions() );
  moncode2.run();
 
  Localisation<2,10,3> moncode3( argc, argv, makeAbout(), makeOptions() );

  //moncode1.run();
  moncode3.run();


}



